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A thin dielectric plate can support conventional slab modes such as leaky Fabry-Perot, guided waves and 
radiation modes. When the plate is coated on a metallic subwavelength disk array, it can result in new 
confined modes, hereby named 'disks scattering slab modes' (DSSMs). By use of a terahertz time-domain 
spectroscopy system, we experimentally verify the existence of two types of DSSMs in terahertz range: one 
we refer to as highly confined disks scattering Fabry-Perot-like modes (DSFPMs) and the other as, weakly 
confined disks scattering radiation-like modes (DSRMs). Spectral characteristics of these confined modes 
are measured and numerically simulated. Particularly, based on the experimental results and a 
dipole-scattering model developed here, we show that the features of the DSFPMs can be tuned by changing 
the filling fraction of the disk array. We believe these results can make important contributions to the 
designs of new terahertz devices, including tunable absorbers and filters. 

I t is well known that surface plasmon polariton-like bound modes, namely spoof surface plasmons (SSPs), can 
I be sustained by periodically structured metal surfaces for longer wavelengths 13 . The most exciting character- 
I istic of SSPs is that their resonant frequencies and confinement can be tuned by the geometric parameters of 
the metallic periodic structures. This feature has recently sparked a wealth of experimental and theoretical works, 
which are aiming at developing new applications, including plasmonic chips 4 , filters 5 7 , and highly sensitive 
imaging 8 . However, SSPs supported by freestanding plasmonic metal structures are weakly confined for most 
frequencies, thus their performance in some of these practical applications is unsatisfactory. In order to resolve 
this problem, a class of highly confined modes was proposed in infrared region, namely hybrid spoof surface 
plasmons (HSSPs) 9 . These modes are supported by a metal-dielectric heterostructure (MDHS) comprised of a 
superstrate of dielectric layer and an ultra thin perforated metal film. The two types of HSSPs are the hybridization 
between SSPs and Fabry-Perot (FP) modes, namely leaky SSP-FP and the hybridization between guided waves 
(GW) and SSPs, referred to as SSP-GW. Based on the concept of HSSPs, practical applications such as tunable 
absorbers and multispectral infrared imaging can be anticipated. 

Unlike the heterostructures mentioned above, we fabricate a set of samples in which hole arrays are switched 
for complementary geometries consisting of metallic disks that are printed on one surface of a thin dielectric slab. 
The reflectance of these samples is measured by a terahertz time-domain spectroscopy (THz-TDS) system (details 
about the system are given in the Methods section) whereby an incident THz beam is normally introduced onto 
the non-arrayed side of the sample. The measured reflection spectra reveals several reflectance dips, that also 
correspond to absorption maxima, just as in the case of heterostructures with hole arrays 9 . However, these 
reflection dips observed here cannot be supported by HSSPs anymore because conventional SSPs do not exist 
in isolated metallic structures 1,2-10 , such as disk/rectangle arrays. Moreover, the observed phenomenon differs 
from both reflection behaviors of metallic disk arrays without a dielectric layer 11 and metallic disk arrays with a 
dielectric layer that has a THz beam illuminating the disk arrayed surface 12 . In both of the above reflectance peaks, 
instead of dips, are usually expected or observed. Moreover, these works show that the reflection resonances are 
attributed to quasi bound states 11 (QBS) inherent to freestanding disk arrays, or dipolar localized surface plas- 
mons 12 (DLSPs) sustained by supported disk arrays. However, for our samples, it is clear that the reflection dips do 
not result from QBS, DLSPs, and HSSPs. The question then remains, how does this phenomenon arise? 

In order to answer this question, we suppose a new class of confined modes, which we shall call disks scattering 
slab modes (DSSMs) originating from the interaction among the scattered fields, diffractions, and conventional 
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slab modes including FP, GW or radiation modes. By analyzing the 
measured results, we find two types of DSSMs. The first exists in 
lower frequencies and its resonances present a redshift behavior with 
increasing filling fraction of the disk arrays. The second type appears 
in higher frequencies and its resonances are nearly independent of 
the geometry of the disk arrays. In order to reveal the physical process 
of these new modes, we perform a numerical simulation by means of 
the rigorous coupled wave analysis (RCWA) method (details about 
the software based on RCWA are given in the Methods section). The 
simulated results indicate that the first type is highly confined, closely 
relating to the FP mode, and hence we name it, disks scattering 
FP-like mode (DSFPM), while the second type is weakly confined, 
relating to the radiation mode, and thus hereby referred to as disks 
scattering radiation-like mode (DSRM). 

From the viewpoint of application, DSFPMs appear more inter- 
esting and valuable due to the ability to be controlled by geometrical 
design of the meta- surface. Hence closer attention is paid to DSFPMs 
than DSRMs in this work. The DSFPM mechanism is explained by 
developing a dipole-scattering model, which shows that this kind of 
resonance results from the coupling between the conventional leaky 
FP modes inherent to low-index dielectric layer and the multiple 
reflections of the ± 1 diffracted beams produced by wave scattering 
on the disk array. All results encountered here may find important 
applications in controlling THz and microwave in both spatial and 
spectral domains. 

Results 

Figure la illustrates a heterostructure comprised of a dielectric layer 
and a meta-surface such as an array of metallic circular disks, where 
the thickness of the dielectric layer, the diameter, and period of the 
disks are denoted by h 0 , a, and p, respectively. We suppose that a 
plane THz wave illuminates the structure at normal incidence, and its 
polarization is parallel to the x axis. Compared to the incident wave- 
length, both p and a are small. In the experiment, we fabricate the 
heterostructure using printed circuit boards (PCB) techniques, as 
shown in Fig. 1. By fixing the reflection spectrum of a silver mirror 
as the reference, the reflectance of the samples is obtained by a THz- 
TDS system, as presented in Fig. lb. The lowest available frequency 
of our THz-TDS for reflection mode is about 0.06 THz. 

We first investigate the effects of the FR4 thickness {h D ) on reflec- 
tance, as illustrated in Fig. 2a. One can see from the figures that the 
existence of metal structures has a large impact on the reflection 
spectra when h 0 is smaller than 770 um. The reflection dips and 
their shift are due to the presence of DSSMs, which is discussed more 
clearly later. We can also deduce from Fig. 2a that the interaction 
among the scattered fields, diffractions, and conventional slab modes 
decreases with increasing h D , derived from the increasing loss of the 
FR4 layer. As a result, when the layer is thicker than 770 urn, 
the excitation of DSSM is weak enough to be neglected. Accordingly, 
the case for h u = 240 um and h D = 500 um will be discussed later. 

In addition to the reflection spectra, the corresponding absorption 
spectra are also shown as reference in Fig. 2b, e, and f. These absor- 
bance are calculated using directly measured reflectance and trans- 
mittance. It should be noted that the lowest credible frequency of our 
THz-TDS for transmission mode is about 0.1 THz, which is higher 
than that of the reflection mode. Thus calculated absorption spectra 
are not accurate below 0.1 THz (the shaded areas in Fig. 2). However, 
one can still see from Fig. 2 that the positions of the reflection dips 
approximately correspond to those of the absorption peaks. As men- 
tioned above, these dips and peaks are associated with DSSMs. 

Next, we measure the reflection spectra of the samples with vari- 
able disk period p and diameter a, respectively (not displayed here). 
We find that the reflectance dips at DSSM's frequencies do not shift 
monotonically with increasing/) or a. This behavior is very different 
from that of metallic disk arrays without any dielectric layer 11 or 
those with thick dielectric layer 12,13 . Compared with these previous 
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Figure 1 | The experiment, (a) Schematic of the metal-dielectric 
heterostructure. Note that Ei = e 3 = 1 (air), e 2 = £ fr4- (b) Schematic of the 
experimental 'set-up'. Ml, M2, and M3 are reflection mirrors. PCA is a 
photoconductive antenna. TA is a THz absorber. OAP1 and OAP2 are off- 
axis parabolic mirrors. Pump beam and probe beam are both femtosecond 
lasers of 800 nm with a time delay between them. The probe beam and THz 
beam are made collinear using ITO glass (transmits optical and reflects 
THz). HRS is high resistance silicon (Si). TD is a THz detector, which 
consists of an EO crystal, a quarter-wave plate, a Wollaston prism, and a 
pair of photodiodes. The inset to Fig. lb is a microscope image of a sample 
corresponding to the heterostructures (bottom view). The samples are 
fabricated by using well-developed printed circuit boards (PCB) 
techniques of area about 2 cm X 2 cm with a copper (Cu) disk array 
printed onto an epoxy glass-cloth laminated layer (FR4). The metallic 
structure is ultrathin (about 17 um) compared to incident wavelength, 
thus its thickness can be neglected. In addition, the dielectric permittivity 
of FR4 is obtained as £ FR4 = 4.4 + 0.35i through its transmittance 
measured by THz-TDS. 

works, there are two important differences in our experiment: First, 
the disk array is covered with a dielectric layer, whose thickness is 
thin but cannot be neglected. Second, the THz wave is incident from 
the side of the air-FR4 interface instead of the air-disk interface. 
Accordingly, the reflection resonances do not lead to spectral peaks 
as in previous works 1012 , however they do produce some reflection 
dips, as illustrated in Fig. 2a, c, and d. Furthermore, the resonant 
frequencies do not shift monotonically with increasing p or a. 

To clearly describe the influence of the geometric parameters 
of the disk array, filling fraction (F = na 2 /4p 2 ) is introduced. 
The THz reflectance or absorbance as a function of F is shown in 
Fig. 2c to f, where the cell gap (p-a) is fixed to 225 pm. As mentioned 
above, there are two types of DSSMs, which lead to reflectance dips 
or absorbance peaks. The spectral resonances denoted by square 
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Figure 2 | The experimental results, (a, b) The effects of FR4 layer thickness on the reflection and absorption spectra. The reflection spectra of 
heterostructure samples (black solid curves) with h D of 240 um, 500 um, 770 um, and 970 urn are measured and compared with the spectra of 
corresponding bare FR4 (red dashed curves), respectively. The geometric parameters of the samples measured here are fixed to p = 800 um, and a = 
500 um. Similar behaviour is observed for different p and a. (c, d) The filling fraction dependence of the reflection spectra of samples with h D = 240 um 
(c) and hu = 500 urn (d), when cell gap (p-a) is fixed to 225 um. (e, f) The absorption spectra corresponding to (c) and (d) respectively. Similar behaviour 
can be obtained for different p-a, which is not displayed here. Note that the calculated absorption spectra are not accurate below 0.1 THz. 



brackets are attributed to resonant coupling to DSFPMs. One can 
observe that these resonances are redshifted when F is increased. It 
means that these modes are sensitive to the geometry of the disk 
array, similar to the leaky SSP-FP resonances 9 . Besides, when F is 
larger than 22.2%, some additional resonances at higher frequencies 
(reflection dips or absorption peaks) can be observed at about 
0.44 THz (Fig. 2c) and 0.38 THz (Fig. 2d), respectively. Similarly 
as in the case of SSP-GW, these resonances are not sensitive to the 
geometry of the metallic structure array, and we attribute them to the 
result of DSRMs, whose field distribution is discussed below. 

Discussion 

By using the RCWA method, the resonant frequencies and field 
distributions of the samples can be obtained. As a specific example, 
we present the transverse distributions of THz fields (\H y \ and |£ z |) in 
one of the samples at the resonant frequencies of 0.19 and 0.43 THz, 
shown in Fig. 3. The two frequencies calculated through RCWA are 
close to the data measured in the experiment. We could infer from 



Fig. 3a and b that the incident THz wave can excite a kind of highly 
confined mode in the MDHS at lower frequencies, such as 0.19 THz 
here. It is clear that the fields of this confined mode are strongly 
concentrated near the disk array, thus its spectral position is strongly 
influenced by the geometric parameters of the array. As mentioned 
in the experimental results, this type of highly confined mode is 
referred to as DSFPM. In addition to DSFPMs, there is another kind 
of weakly confined mode labeled DSRM, resulting from the inter- 
action between the fields scattered by disks and the conventional 
radiation mode supported by the FR4 waveguide layer. The field 
distribution of DSRMs is mainly concentrated near the air-FR4 inter- 
face, shown in Fig. 3c. Therefore, the spectral positions of DSRMs are 
not sensitive to the geometry of the disk arrays. The numerical ana- 
lysis above supports our experimental results very closely. From the 
viewpoint of application, the DSFPMs are more interesting and valu- 
able due to the ability to be controlled by the geometrical design of 
the meta-surface. Thereby the DSRMs will not be discussed in detail 
below. 



SCIENTIFIC REPORTS | 3 : 2493 | DOI: 10.1038/srep02493 



3 



200. 



1001 



a \H y \ /=0.19 THz 



-100 



-200 




-300 -200 -100 0 100 200 300 

a- (urn) 



|1.4 


200, 


-1.2 


100 


1 

0.8 


0 

1-100 


0.6 


-200 


0.4 




0.2 


-300 
-400 



b 






\Ez 


/=0.19 THz 



-300 -200 -100 0 100 200 300 

*(um) 



3 

2.5 
2 

1.5 
1 

0.5 



200 








200r 




-1.2 


100 




1 








0 




0.8 
0.6 


1-100 




0.4 


-200 




0.2 


-300 
-400 



-300 -200 -100 0 100 200 300 

x (urn) 




-300 -200 -100 0 100 200 300 

x (urn) 




slice 




F=36.4% 



Figure 3 | Near field distribution, (a, b) The transverse distribution of the magnetic field ( \H y \ ) and electric field (lEjl) for the case of/= 0.19 THz, 
which corresponds to the DSFPM. (c, d) The transverse distribution of the magnetic field ( I H y I ) and electric field ( IEJ ) for the case of / = 0.43 THz, 
which corresponds to the DSRM. The bottom panels display the unit cell and slice of the MDHS calculated in RCWA software. In the numerical 
computation, the geometric parameters are chosen as F = 36.4%, h D = 240 fim, and p-a = 225 um. 



To explain the red-shift behavior of the DSFPMs, a laconic for- 
mulation, not requiring special function evaluations, has been estab- 
lished by considering the MDHS as an asymmetric Fabry-Perot 
cavity and applying the dipole-scattering approach 14 . Under normal 
incidence, the total reflection matrix of the structure shown in Fig. la 
can be expressed as: 



•"total — r 12 +t2lPr p 23ZPt2i, 



Z=fI-Pr 21 Pr, 



2irr p23 j 



(la) 
(lb) 



where P describes plane-wave propagation across the FR4 layer, 
(ty) is the reflection (transmission) matrix at the rth-jth medium 
interface, and r p23 represents the reflection matrix of the periodic 
disk array. For a freestanding disk array, when the incident wave- 
length is commensurate with the period of the disks p (or multiple 
period mp) scattering of the incident wave by the disk array produces 
diffracted beams that are evanescent above their diffraction thresh- 
old. However, in the presence of a coated layer these evanescent 
diffracted beams may undergo multiple reflections in the asymmetric 
cavity formed between the air-FR4 interface and the array, thereby 
contributing to the confinement of THz wave. It should be noted that 



since the external field is polarized along the x axis, the diffracted 
beams parallel to x are not expected to play a significant role 15 . 
Therefore, here we include some diffraction beams with parallel wave 
vectors g m = (mgjy, in which = 2n/p and m = 0, ±1, ±2, 
±m MAX . Accordingly, r,j can be expressed as a (2m MAX + 1) X 
(2m MAX + 1) matrix with coefficients r,y (mm ', relating the 
incident beam m' to reflected beams m. By using Fresnel's formula 16 , 
we find 



(mm') _ c 
'12 u mm 



-eik 2 /e 2 ki 
-s 1 k 2 /e 2 ki ' 



(mm') 



1 e 2^1zmAl^2zm 
1 +£ 2 klzm/ 2l 



and the transmission matrix satisfies t^ n 



(2a) 

(2b) 
■'». In 



Eq. (2), k km = ijkj-imgi) 2 (m = 0, ±1, ±2, ±m MAX ) is the z 

component of the wave vector (fc; = k 0 nj), where A: 0 is the free-space 
wave vector, and n ( is the refractive-index of medium Z. Moreover, the 
coefficients of matrix P can be expressed as p* mm '' = <5 mm - 
exp(ik 2zm h D ), respectively. 
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Next, applying the dipole-scattering approach, each disk can be 
assumed to respond with an induced dipole, thus the reflection coef- 
ficients of r p2 3 are obtained as 



(mm') _ 

'p23 " 



2nik\ 



P 2 k 2 



-). 



(3) 



where a E is the electric polarizability of effective dipoles, and G E is the 
dipole-dipole interaction dyadic. Due to the presence of the coated 
dielectric layer, the electric polarizability a E is different from that of 
freestanding disk array. For metals in low frequency regime, perfect 
electric conductor (PEC) approximation leads to simple analytical 
expression of the polarizability, 

* E a 3 £2 1 3 ' (4j 

in which t, is deduced as (£2-£3)/(£2+ £ 3)> which describes the asym- 
metric surrounding environment of the disk array. Particularly, C = 
0 means that the array of disks is placed in a homogeneous medium 
without coated layer or substrate. For the case of ). 2 >p and a 2 ~ p, 
another important parameter G E reduces to 



47i 2 v / 2 



P 



+K- 



2nk 2 



(5) 



where C(l 2 /p) is a smooth fit to the remaining nondivergent terms 11 . 
Then, based on the equations above, the zero-order coefficient Z <00) of 
the multipass matrix Z can be easily obtained. Finally, upon inspec- 
tion of the denominator of Z (m , we can express the phase matching 
condition of the DSFPMs in a simple form, 



(6a) 



S m = r 2 " m) r ( ^ ) exp(2ik 2zm h D ), 



(6b) 



where S m describes multiple reflections of the diffraction beam m, q is 
an positive integer number, and/ R is the resonant frequency of the 
DSFPM q. To investigate the effect of m MAX on the accuracy of our 
model, the calculation results for m MAX = 0, 1, and 2 have been 
shown in Fig. 4a, in which the experimental results are also given. 
In fact, if only the specular beam is considered in our model (m = 
m MAX = 0), Eq. (6) can be written as arg{S 0 } = 2qn, which is the 
conventional Fabry-Perot formula. Therefore it is inaccurate to neg- 
lect higher order diffraction beams, as presented in Fig. 4a. One can 
see from the calculation results that we can get a minimal analytical 
formulation with enough accuracy by fixing m MAX = 1, which indi- 
cates that three diffraction beams with parallel wave vectors g 0 = 0 
and g±i = ± {2nlp)y are included in our theory. In this case, the 
phase matching condition of the DSFPMs can be expressed as 



arg{S_ i + S 0 + Si } /=fe = 2qn 



S 0 = r! 



(oo) ,(oo) 

21 



P 23 exp(2* 2 /z D ), 



Si =S_! = r 2 \ 1) r ( p 1 2 °^exp(2ik 2z ih D ). 



(7a) 
(7b) 
(7c) 



Based on the discussion above, Eq. (7) gives us an insight into the 
physical mechanism of DSFPMs. The multiple reflections of the 
specular beam (S 0 ) lead to a conventional FP behaviour in the het- 
erostructure. Meanwhile, this process can receive important contri- 
butions from the multiple reflections of diffracted beams (S ±1 ) with 
g ±1 = ± (27t/p)y via momentum exchange with the disk lattice. 
Compared to the conventional FP formula arg{S 0 } = 2qn, Eq. (7) 
gives a FP-like phase matching condition, thus the highly confined 
DSSMs discussed here are referred to as disks scattering FP-like 
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Figure 4 | The theoretical discussion, (a, b) The filling fraction dependence of the DSFPM frequencies for samples shown in Fig. la. Note that q = 1 
corresponds to dip 1 in Fig. 2c and d, while q = 2 represents dip 2 in them, (c, d) The results of metal-dielectric square patch arrays, for which the filling 
fraction is defined as F = a 2 /p 2 . When the disk array on the bottom of the sample is replaced by square patch array, the real part of (a £ ) _1 should be 
changed, and an analytical expression similar to Eq. 4 is unavailable. However, we can still estimate that Re{(a; £ ) _1 ( oc a~ 3 (l + £)/e 2 . The insets to Fig. 4 are 
microscope images of unit cells of the samples (bottom view) . In all figures, the theoretical results are represented by solid lines, and the experimental 
results are denoted by dashes and symbols. Note that the cell gap (p-a) is fixed to 225 um in these discussion. 
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modes (DSFPMs). It should be noted that there are also contributions 
from higher order diffracted beams due to scattering in this physical 
process, but these contributions are insignificant and thus can be 
neglected. Furthermore, if the disk arrays in the MDHS are replaced 
with hole arrays, we could extend our theory to the case of HSSPs by 
applying Babinet's principle 14 to r p23 . Then the HSSP phase matching 
condition obtained in Ref. 9 can be replaced by a similar form. 

Through Eq. (7), we calculate the resonant frequencies of the 
DSFPMs with increasing F for the case of h D = 240 urn and 
500 urn, respectively. These theoretical results are compared to the 
experimental data as shown in Fig. 4a and b. It is clear that the 
changing tendency of the theoretical prediction is in accordance with 
the measurement. Particularly, for h v = 240 um, experimental data 
and theory are closely matched. Meanwhile, we also find that the 
prediction of our model is not accurate for relatively larger filling 
fractions. This is because the dipolar scattering theory used to obtain 
r p2 3 does not work well when the size of the array is far away from the 
subwavelength limit. Moreover, if multipole approximation and 
complete diffraction beams 915 are incorporated into the model, the 
accuracy of prediction can be improved, but a particularly compli- 
cated formulation will be encountered. 

Furthermore, we can deduce from Fig. 4a and b that the DSFPMs 
of q = 1 is the dominant factor in the resonance-shift behavior for h v 
= 240 urn, while the modes of q = 1 and 2 are predominant for the 
case of h D = 500 um. Upon inspection of equation (7), one finds that 
the phase matching condition is determined not only by r p23 but also 
by r 2 i and P 2 , thus the effect of {II a E -Ge)" 1 on r tota i decreases, which 
is different from the case of freestanding disk arrays (r, otal = r p23 ). 
Therefore, the resonant frequencies do not shift monotonically with 
increasing a or p as previous works 6,1113 , and the introduction of F is 
suitable here to describe the additional phase caused by the array- 
dielectric interface. Finally, we extend the investigation to the case of 
metal- dielectric square patch array, as shown in Fig. 4c and d, which 
support our previous analysis well. Based on the discussion above, 
the experimental results are explained, and the existence of DSSMs is 
demonstrated. 

It is well known that lattice modes can exist in metallic nanodisk 
arrays and will hybridize with the localized surface plasmons (LSPs) 
sustained by each of the metal disks 17 20 . However, such hybridized 
modes are totally different from DSFPMs, because in the THz regime 
real plasmons, including LSPs, do not exist in the metal disks, which 
are usually treated as PEC disks for THz wave. Meanwhile, the pres- 
ence of the FR4 layer can also suppress the so called lattice modes 
supported by the disk arrays. In addition, it should be noted that a 
new class of hybridized confined modes, i.e. LSP-FP modes, sus- 
tained by an isolated nanoparticle-cavity system, have been reported 
in the visible range recently 21 . Although both of LSP-FP modes and 
DSFPMs are highly confined FP-like modes, LSP-FP modes can 
either enhance or suppress the scattering intensity slightly while 
DSFPMs can strongly suppress the reflectance. 

In summary, we have experimentally and theoretically demon- 
strated that a heterostructure, comprised of a dielectric layer and a 
metallic subwavelength disk array, can support a new class of con- 
fined modes, namely disks scattering slab modes, in THz domain. 
These modes originate from the interaction among the scattered 
fields, diffractions, and conventional slab modes such as FP or radi- 
ation modes. Thus they include two types: DSFPMs and DSRMs. 
Particularly, the DSFPMs are highly sensitive to the geometry of 
the disk arrays, and their resonances are redshifted when the array's 
filling fraction is increased. Based on these results, new THz applica- 
tions in spatial and spectral control can be expected, such as tunable 
absorbers and filters. 

Methods 

The experimental setup used in this work is based on a commercial THz-TDS 
system, which is known as Z3. It measures far-infrared spectra between 0.06 THz to 



3.5 THz {reflection mode). The peak dynamic range is greater than 70 dB. The 
system is pumped using a pulsed fiber laser with a central frequency at 800 nm, 
average power of 100 mW and a pulse width of 100 fs. Oscillators with up to 1 W 
average power can improve system performance by as much as 20 dB. The laser 
beam is split into a pump and probe beam used for THz wave generation and 
detection. A photoconductive antenna (fabricated on a LT-GaAs wafer) is used to 
generate THz waves. The time delay between the pump and probe is controlled by 
dual retro-reflectors driven by a servo. In the THz detector, the probe beam travels 
collinearly with the THz beam in the EO crystal, where the probe beam is modu- 
lated by the electric field of the THz radiation via an electro-optic effect. A quarter- 
wave plate, a Wollaston prism and a pair of photodiodes are assembled for the 
balanced detection of the probe beam. 

The experimental samples are designed through Protel 99SE, and fabricated by 
PCB techniques. We monitor and measure the geometric parameters of the samples 
by using of a metaloscope. Before the THz reflection experiment, each of the samples 
is cleaned by acetone. From the reflected/transmitted waveforms measured directly 
by the THz-TDS, we obtain the reflection/transmission spectra of all samples by using 
Fourier transformation. All measurements are performed with dry air purging (0.5% 
humidity). The absorbance is calculated through these measured results, and then 
corrected by removing the baselines. Note that the lowest available frequency of our 
THz-TDS for transmission mode is about 0.1 THz, which is different from that of 
reflection mode, thus the calculated absorption spectra are not credible below 
0.1 THz. 

To obtain the resonant frequencies and near field distributions of the MDHS 
samples, the commercially available RCWA software DiffractMOD is used. In the 
simulation, the permittivity of the copper disks is fixed as c M - -3 X 10 4 + 10 5 i 22 , 
while that of the FR4 layer is chosen as e FR4 — 4.4 + 0.35i. Periodic boundary 
conditions are imposed in x and y directions to consider periodic arrangement of the 
unit cells shown in the bottom panels of Fig. 3. 
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